#This script includes the code necessary to replicate the Tables and Figures from Benn, "Presidential Partisanship and Regulatory Review"

pacman::p_load(tidyverse, stargazer, plm, lfe, bife, foreign, sandwich, lmtest)
load("./Replication.Rdata")

#Table 1 ----
# Baseline specifications, asymmetric partisanship hypothesis

#Split data by partisanship
regdem <- regdata[regdata$Republican == 0, ]
regrep <- regdata[regdata$Republican == 1, ]

agdem <- merged[merged$Republican == 0, ]
agrep <- merged[merged$Republican == 1, ]

#Baseline specifications
lm_b1 <- felm(data = regdata, formula = Change ~ Republican + RCL + party_interact_RCL + Attention + review.length| 0 | 0 | AGENCY_CODE)
lm_b2 <- felm(data = merged, formula = Change ~ Republican + RCL + party_interact_RCL + Attention + review.length + words + significance| 0 | 0 | AGENCY_CODE)
lm_b3 <- felm(data = regdem, formula = Change ~ RCL + Attention + review.length| 0 | 0 | AGENCY_CODE)
lm_b4 <- felm(data = regrep, formula = Change ~  RCL + Attention + review.length| 0 | 0 | AGENCY_CODE)
lm_b5 <- felm(data = agdem, formula = Change ~ RCL + Attention + review.length + words + significance| 0 | 0 | AGENCY_CODE)
lm_b6 <- felm(data = agrep, formula = Change ~ RCL + Attention + review.length + words + significance| 0 | 0 | AGENCY_CODE)

#Table output
stargazer(lm_b1, lm_b2, lm_b3, lm_b4, lm_b5, lm_b6,
          header = FALSE, 
          type='latex', 
          title = "Ideological Alignment Hypothesis",
          covariate.labels = c("Republican", "Agency Ideology (RCL)", "Interaction", "Attention", "Review Length", "Complexity", "Significance", "Constant"),
          model.names = FALSE,
          column.sep.width = "1pt",
          add.lines = list(c("Data", "Regulatory", "Agenda", "Reg. (Dem.)", "Reg. (Rep.)", "Agenda (Dem.)", "Agenda (Rep.)")),
          omit.stat = c("ser", "rsq"))


# Table 2 ----
# Liberal agencies hypothesis

lm_l1 <- felm(data = regdata, formula = Change ~  RCL + Republican + Attention + review.length | 0 | 0 | AGENCY_CODE)
lm_l2 <- felm(data = merged, formula = Change ~ RCL + Republican + Attention + review.length + words + significance| 0 | 0 | AGENCY_CODE)
lm_l3 <- felm(data = regdata, formula = Change ~  RCL + Attention + review.length | President | 0 | AGENCY_CODE)
lm_l4 <- felm(data = merged, formula = Change ~ RCL + Attention + review.length + words + significance| President | 0 | AGENCY_CODE)
lg1 <- glm(data = regdata, formula = Change ~  RCL + Republican + Attention + review.length, family = "binomial")
lg2 <- glm(data = merged, formula = Change ~ RCL + Republican + Attention + review.length + words + significance, family = "binomial")

stargazer(lm_l1, lm_l2, lm_l3, lm_l4, lg1, lg2,
          title = "Liberal Agencies Hypothesis",
          font.size = "small",
          omit.stat = c("ser", "ll", "aic", "rsq"),
          header = FALSE,
          model.names = FALSE,
          covariate.labels = c("Agency Ideology (RCL)", "Republican", "Attention", "Review Length", "Complexity", "Significance", "Constant"),
          add.lines = list(c("Model", "OLS", "OLS", "OLS", "OLS", "Logit", "Logit"),
                           c("President FE", "N", "N", "Y", "Y", "N", "N"),
                           c("Data", "Regulatory", "Agenda", "Regulatory", "Agenda", "Regulatory", "Agenda")))


# Table 3 ----
m1 <- felm(data = regdata, formula = Change ~  Republican + Attention + review.length | AGENCY_CODE | 0 | 0)
m2 <- felm(data = merged, formula = Change ~  Republican + Attention + review.length + words + significance | AGENCY_CODE | 0 | 0)
m3 <- felm(data = regdata, formula = Change ~  Republican + RCL + dim1est + dim2est + Attention + review.length| 0 | 0 | 0)
m4 <- felm(data = merged, formula = Change ~  Republican + RCL + dim1est + dim2est + Attention + review.length + words + significance | 0 | 0 | 0)
m5 <- glm(data = regdata, formula = Change ~  Republican + RCL + dim1est + dim2est + Attention + review.length, family = "binomial")
m6 <- glm(data = merged, formula = Change ~  Republican + RCL + dim1est + dim2est + Attention + review.length + words + significance, family = "binomial")

r1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC1"))
r2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC1"))
r3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC1"))
r4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC1"))
r5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC1"))
r6 <- coeftest(m6, vcov = vcovHC(m6, type = "HC1"))

stargazer(m1, m2, m3, m4, m5, m6,
          se = list(r1[,2], r2[,2], r3[, 2], r4[, 2], r5[, 2], r6[, 2]),
          title = "Republican Presidents Hypothesis",
          font.size = "small",
          omit.stat = c("ser", "ll", "aic", "rsq"),
          model.names = FALSE,
          header = FALSE,
          covariate.labels = c("Republican", "Agency Ideology (RCL)", "Independence DM", "Independence PD", "Attention", "Review Length", "Complexity", "Significance", "Constant"),
          add.lines = list(c("Model", "OLS", "OLS", "OLS", "OLS", "Logit", "Logit"),
                           c("Data", "Regulatory", "Agenda", "Regulatory", "Agenda", "Regulatory", "Agenda"),
                           c("Agency FE", "Y", "Y", "N", "N", "N", "N")))

#Table 4 ----
lg_b1 <- glm(data = regdata, formula = Change ~ Republican + RCL + party_interact_RCL + Attention + review.length, family = "binomial")
lg_b2 <- glm(data = merged, formula = Change ~ Republican + RCL + party_interact_RCL + Attention + review.length + words + significance, family = "binomial")
lg_b3 <- glm(data = regdem, formula = Change ~ RCL + Attention + review.length, family = "binomial")
lg_b4 <- glm(data = regrep, formula = Change ~  RCL + Attention + review.length, family = "binomial")
lg_b5 <- glm(data = agdem, formula = Change ~ RCL + Attention + review.length + words + significance, family = "binomial")
lg_b6 <- glm(data = agrep, formula = Change ~ RCL + Attention + review.length + words + significance, family = "binomial")

stargazer(lg_b1, lg_b2, lg_b3, lg_b4, lg_b5, lg_b6,
          header = FALSE, 
          title = "Ideological Alignment Hypothesis, Logit Specifications",
          covariate.labels = c("Republican", "Agency Ideology (RCL)", "Interaction", "Attention", "Review Length", "Complexity", "Significance", "Constant"),
          model.names = FALSE,
          column.sep.width = "1pt",
          add.lines = list(c("Data", "Regulatory", "Agenda", "Reg. (Dem.)", "Reg. (Rep.)", "Agenda (Dem.)", "Agenda (Rep.)")),
          omit.stat = c("ll", "aic"))

#Table 5 ----
m1 <- felm(data = regdata, formula = Change ~  Republican + Attention + review.length| AGENCY_CODE | 0 | AGENCY_CODE)
m2 <- felm(data = merged, formula = Change ~  Republican + Attention + review.length + words + significance | AGENCY_CODE | 0 | AGENCY_CODE)
m3 <- felm(data = regdata, formula = Change ~  Republican + RCL + dim1est + dim2est + Attention + review.length| 0 | 0 | AGENCY_CODE)
m4 <- felm(data = merged, formula = Change ~  Republican + RCL + dim1est + dim2est + Attention + review.length + words + significance | 0 | 0 | AGENCY_CODE)

stargazer(m1, m2, m3, m4,
          title = "Republican Presidents Hypothesis, Clustered SEs",
          font.size = "small",
          header = FALSE,
          omit.stat = c("rsq", "ser"),
          model.names = FALSE,
          covariate.labels = c("Republican", "Agency Ideology (RCL)", "Independence DM", "Independence PD", "Attention", "Review Length", "Complexity", "Significance", "Constant"),
          add.lines = list(c("Model", "OLS", "OLS", "OLS", "OLS"),
                           c("Data", "Regulatory", "Agenda", "Regulatory", "Agenda"),
                           c("Agency FE", "Y", "Y", "N", "N")))

#Table 6 ----
#Table 1 columns 1 and 2
a1 <- felm(data = regdata, formula = Change ~ Republican + RCL + party_interact_RCL| 0 | 0 | AGENCY_CODE)
a2 <- felm(data = merged, formula = Change ~ Republican + RCL + party_interact_RCL + words + significance| 0 | 0 | AGENCY_CODE)
a3 <- felm(data = regdata, formula = Change ~  RCL + Republican | 0 | 0 | AGENCY_CODE)
a4 <- felm(data = merged, formula = Change ~ RCL + Republican + words + significance| 0 | 0 | AGENCY_CODE)
a5 <- felm(data = regdata, formula = Change ~  Republican + RCL + dim1est + dim2est| 0 | 0 | 0)
a6 <- felm(data = merged, formula = Change ~  Republican + RCL + dim1est + dim2est + words + significance | 0 | 0 | 0)

stargazer(a1, a2, a3, a4, a5, a6,
          title = "Robustness check omitting Attention and Review Length",
          header = FALSE,
          font.size = "small",
          omit.stat = c("rsq", "ser"),
          covariate.labels = c("Republican",
                               "Agency Ideology (RCL)",
                               "Interaction",
                               "Independence DM",
                               "Independence PD",
                               "Complexity",
                               "Significance",
                               "Constant"),
          add.lines = list(c("Replication", "Table 1", "Table 1", "Table 2", "Table 2", "Table 3", "Table 3"),
                           c("Data", "Regulatory", "Agenda", "Regulatory", "Agenda", "Regulatory", "Agenda"))
)

#Table 7 ----

e1 <- felm(data = regdata, formula = Change ~ Republican + RCL + party_interact_RCL + Attention + review.length + EO| 0 | 0 | AGENCY_CODE)
e2 <- felm(data = regdata, formula = Change ~  RCL + Republican + Attention + review.length + EO | 0 | 0 | AGENCY_CODE)
e3 <- lm(data = regdata, formula = Change ~  Republican + RCL + dim1est + dim2est + Attention + review.length + EO)
e4 <- glm(data = regdata, formula = Change ~  Republican + RCL + dim1est + dim2est + Attention + review.length + EO, family = "binomial")

er3 <- coeftest(e3, vcov = vcovHC(e3, type = "HC1"))
er4 <- coeftest(e4, vcov = vcovHC(e4, type = "HC1"))

stargazer(e1, e2, e3, e4,
          title = "Robustness check controlling for EO 12866 regime",
          header = FALSE,
          font.size = "small",
          model.names = FALSE,
          omit.stat = c("rsq", "ser", "ll", "aic", "f"),
          se = list(e1$cse, e2$cse, er3[,2], er4[,2]),
          covariate.labels = c("Republican",
                               "Agency Ideology (RCL)",
                               "Interaction",
                               "Independence DM",
                               "Independence PD",
                               "Attention",
                               "Review Length",
                               "EO 12866",
                               "Constant"),
          add.lines = list(c("Replication", "Table 1", "Table 2", "Table 3", "Table 3"),
                           c("Model", "OLS", "OLS", "OLS", "Logit"))
)

# Figure 1 ----
#A function to input president and return coefficient on Republican with lower and upper CIs (full data)

getcoef_full <-
    function(president){
        model <- felm(data = regdata[regdata$President != president,], formula = Change ~ Republican | AGENCY_CODE | 0 | 0)
        coef <- round(coef(model)[[1]], 2)
        se <- coeftest(model, vcov = vcovHC(model, type = "HC1"))[1, 2]
        lb <- round(coef - 1.96*se, 2)
        ub <- round(coef + 1.96*se, 2)
        return(c(president, coef, lb, ub))
    }


coefs.full <- lapply(X = as.vector(levels(as.factor(regdata$President))), FUN = getcoef_full) %>%
    do.call(rbind.data.frame, .)
names(coefs.full) <- c("President", "Coefficient", "Lower", "Upper")
coefs.full <- 
    coefs.full %>%
    mutate(Lower = as.numeric(Lower), 
           Upper = as.numeric(Upper),
           Coefficient = as.numeric(Coefficient))

#Coefficient plot for full dataset
ggplot(coefs.full, aes(President, Coefficient)) +
    geom_point() +
    geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.1) +
    theme_bw() +
    labs(x = "Omitted President")+
    geom_hline(aes(yintercept = 0, linetype = "longdash")) +
    scale_x_discrete(limits = c("Reagan", "HWBush", "Clinton", "WBush", "Obama", "Trump"),
                     labels = c("Reagan", "H.W. Bush", "Clinton", "Bush", "Obama", "Trump")) +
    guides(linetype = FALSE)


#Table 8 ----
f1 <- felm(data = regdata, formula = Change ~  as.factor(President) | AGENCY_CODE | 0 | 0)
f2 <- felm(data = merged, formula = Change ~  as.factor(President) | AGENCY_CODE | 0 | 0)
f3 <- felm(data = regdata, formula = Change ~ as.factor(President) + Attention + review.length + EO | AGENCY_CODE | 0 | 0)
f4 <- felm(data = merged, formula = Change ~ as.factor(President) + Attention + review.length | AGENCY_CODE | 0 | 0)
f5 <- lm(data = regdata, formula = Change ~ as.factor(President) + RCL + dim1est + dim2est + Attention + review.length)
f6 <- lm(data = merged, formula = Change ~ as.factor(President) + RCL + dim1est + dim2est + Attention + review.length + words + significance)

rr1 <- coeftest(f1, vcov = vcovHC(f1, type = "HC1"))
rr2 <- coeftest(f2, vcov = vcovHC(f2, type = "HC1"))
rr3 <- coeftest(f3, vcov = vcovHC(f3, type = "HC1"))
rr4 <- coeftest(f4, vcov = vcovHC(f4, type = "HC1"))
rr5 <- coeftest(f5, vcov = vcovHC(f5, type = "HC1"))
rr6 <- coeftest(f6, vcov = vcovHC(f6, type = "HC1"))

stargazer(f1, f2, f3, f4, f5, f6,
          se = list(rr1[,2], rr2[,2], rr3[, 2], rr4[, 2], rr5[, 2], rr6[, 2]),
          title = "Presidential Administration Fixed Effects",
          header = FALSE,
          font.size = "small",
          omit.stat = c("rsq", "ser", "f"),
          covariate.labels = c("H. W. Bush",
                               "Obama",
                               "Reagan",
                               "Trump",
                               "W Bush",
                               "Agency Ideology (RCL)",
                               "Independence DM",
                               "Independence PD",
                               "Attention",
                               "Review Length",
                               "EO 12866",
                               "Complexity",
                               "Significance",
                               "Constant"),
          model.names = FALSE,
          add.lines = list(c("Agency FE", "Y", "Y", "Y", "Y", "N", "N"),
                           c("Data", "Regulatory", "Agenda", "Regulatory", "Agenda", "Regulatory", "Agenda")))


